% by ylx,clq on 6/15
% Version 2.0
%% for fig 10 12 14

function [prctages, cost_struct, profits_prct] = figure1x_even(j_all_all)
   
    global I J ci cI r CURRENT s_time S
    ci_origin = ci;
    cI_origin = cI;
    r_origin = r;
    

    cI = 100;
    r = 100;
    cis = 40:10:90;
    N = size(j_all_all, 2); % 可否这样写省去N参数?
    len_seq = size(j_all_all, 1); %认为j_all_all为seq行向量的组合
    len_cis = size(cis, 2);
    index = zeros(1, len_seq); % 每个seq最优利润病人的总数
    cost_struct = zeros(3, I); %收入减去两个成本的矩阵
    opt_profits = zeros(1, N);  %同个seq的最高利润组合
    opt_costruct = zeros(3, N);  %所对应的cost structure
    max_profits = zeros(1, len_seq); %不同seq的最高利润
    max_costructs = zeros(3, len_seq); % 所对应的cost structure
    profits_prct = zeros(3, len_cis); % 每个ci值对应的10,50,90分位的利润
    costructure = zeros(3, len_cis);  % 每个ci值对应的平均利润组成
    i_star = zeros(1, N);
    slot_patientno = zeros(I, len_seq);
    prctages = zeros(I, len_cis);%记录每个ci值的slot分配
    %对每个ci，跑2500次，记录期望利润，利润part，slot
    for jj = 1:len_cis%依次更改ci的值
        % 这里大段代码重复 未来可以打包成函数
        % 代码来自 prop_policy.m
        ci = cis(jj);
        for ii = 1:len_seq % 依次更改j_all
            fprintf('\t fig1x_even -> %d/%d * %d/%d\n', jj, len_cis, ii, len_seq);
            if mod(ii, 10) == 0
                fprintf('\t\t Total time = %.1f\n', ftimer()-s_time);
            end
            j_all = j_all_all(ii, :);
            [best_is, opt_profits, EYI, EYin, opt_costruct,place] = prop_policy2(j_all, N);
            %记录每个j_all利润最高时的利润，利润part，slot分配情况
            max_profits(1,ii) = opt_profits(place);%记录利润
            max_costructs(:,ii)=opt_costruct(:,place);%记录利润组成
            slot_patientno(:,ii)=sum(S,2);%记录slot分配
        end    
        prctages(:,jj) = sum(slot_patientno, 2)/sum(sum(slot_patientno));%记录每个ci值的slot分配
        profits_prct(:, jj) = prctile(max_profits,[10; 50; 90]);%记录每个ci值的利润10，50，90分位
        costructure(:, jj) = mean(max_costructs,2);%记录每个ci的平均利润组成
    end
    
    figure(10) 
    for i = 1:I
        if mod(i,2)==0
            plot(cis, 100*prctages(i, :), '--',"LineWidth",2);
        else
            plot(cis, 100*prctages(i, :),"LineWidth",2);
        end
        hold on
    end
    xlabel("c_i");
    ylabel("Percentage of patients assigned");
    yticks([6:2:18]);
    xticks([40:10:90]);
    legend('1','2','3','4','5','6','7','8','Location','Best')
    ylim([6 18]) 
    grid on
    % <<< 

    % Fig 12
    figure(12)
    for i = 1:3
        plot(cis, profits_prct(i, :),"LineWidth",2);
        hold on
    end
    xlabel("c_i");
    ylabel("Expected Profit");
    xticks([40:10:90]);
    yticks([1000:50:1450]);
    legend("10 prctile","50 prctile" ,"90 prctile","location","Best");
    ylim([1000 1450]);
    grid on
    %input("Press ENTER to continue.\n");


    % Fig 14
    figure(14)
    for i = 1:3
        plot(cis,costructure(i,:),"LineWidth",2);
        hold on
    end
    xlabel("c_i");
    ylabel("Components of expected profit for various c_i");
    xticks([40:10:90]);
    yticks([0:500:2500]);
    legend("r E[\Sigma_iX_i]","c_i E[\Sigma_iY_i]" ,"c_I E[\Sigma_IY_I]","location","Best");
    ylim([0 2500]);
    grid on
    %input("Press ENTER to continue.\n")


    % >>> 恢复global变量的值（如果后续没用到原来的值，不写也行）
    ci = ci_origin;
    cI = cI_origin;
    r = r_origin;
    % <<<
    return;


